Load in libraries

library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(DT)
library(ggpubr)
library(cowplot)
source("scripts/utils/utils.R")
modernacolors()
## <environment: R_GlobalEnv>

Load data

# load ELISpot data
temp <- read_excel("processed_data/ELISpot.xlsx")

dp <- temp[!is.na(temp$counts), ] %>%
  filter(`low cell number (<0.15M)`== "no",
         day != "143")%>%
  mutate(counts = counts + 0.001)

d <-
  dp %>%
  filter(interval != "PBS")%>%
  mutate(animal_id = factor(animal_id),
         interval = factor(interval),
         day = factor(day))

Observed data visualization

dp %>%
  filter(interval != "PBS") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
  ggplot(aes(x = interval, y = counts)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "red",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  facet_wrap(tissue ~ coating, scales = "free") +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "antigen specific \n ASC counts/million cells",
    color = "",
    linetype = "",
    title = "Antigen specific spleen and BM results"
  ) +
  theme(
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "bottom"
  )

Model fitting for spleen data

#-------------------------------------------------------------------------------
# Fit the model for spleen data 
#-------------------------------------------------------------------------------
# boxcox function to fit boxcox regression
boxcox_fit <- function(dat) {
  require(MASS)
  if (length(unique(dat$day)) == 1) {
    out <- boxcox(counts ~ interval, data = dat)
    optimal <- out$x[which.max(out$y)]
    bctran <- make.tran("boxcox", optimal)
    lm_fit <- with(bctran,
                   lm(linkfun(counts) ~ interval, data = dat))
  } else{
    out <- boxcox(counts ~ interval * day, data = dat)
    optimal <- out$x[which.max(out$y)]
    bctran <- make.tran("boxcox", optimal)
    lm_fit <- with(bctran,
                   lm(linkfun(counts) ~ interval * day, data = dat))
  }
  return(list(lm_fit = lm_fit, lambda = optimal))
}

# box cox model fit for spleen data
d_tissue <- subset(d, tissue == "spleen") %>%
  filter(!(coating == "NTD" & interval == "0"))

boxcox_fits <-
  d_tissue %>%
  split(.$coating) %>%
  map( ~ boxcox_fit(.x))

lm_fits <-
  purrr::transpose(boxcox_fits)$lm_fit

# obtain optimal lambda for boxcox transformation
lambdas <-
  purrr::transpose(boxcox_fits)$lambda %>% unlist()

assumption_plots <-
  imap(
    lm_fits ,
    ~ assumption_check(.x, .x$model$`linkfun(counts)`) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0("Spleen ", .y, "-specific model residual diagnostics "),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )

# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
  emmeans(lm_grid, specs = ~ interval, type = "response") %>%
  as_tibble()

Residual diagnostics for spleen data

ggarrange(plotlist = assumption_plots)

#-------------------------------------------------------------------------------
# Figure 4A. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
# fold change figure
foldchange1 <-
  emmt_dat %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(facet = "1 week post-boost") %>%
  ggplot(aes(x = interval , y = response, color = interval)) +
  geom_point(position = position_dodge(width = 0.5), size = 2)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 0.5),
    size = 1
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  #ggtitle("S2P specific ASCs") +
  facet_grid( ~ facet) +
  labs(x = "", y = "Estimated counts/million cells") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(6, 7)],
             color = "darkcyan",
             lty = "dashed") +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

# boxplot
box1 <-
  dp %>%
  filter(interval != "PBS", tissue == "spleen", coating == "S2P") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(facet = "1 week post-boost") %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
  ggplot(aes(x = interval, y = counts, color = interval)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "grey",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  facet_grid( ~ facet) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "Counts/million cells",
    color = "",
    linetype = "",
    title = "S2P specific ASCs"
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
    )
  ) +
  theme(
    axis.text = element_text(size = 12),
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "None",
  )

Supplementary Table 3. Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing intervals.

#-------------------------------------------------------------------------------
# Supplementary Table 3
#-------------------------------------------------------------------------------
emm_elispot <-
  emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
          specs = ~ interval,
          type = "response")

set.seed(100)
sum_elispot <-
  contrast(
    regrid(emm_elispot, transform = "log10"),
    "pairwise",
    infer = TRUE,
    adjust = "mvt",
    ratio = FALSE,
    level = .95
  ) %>% as_tibble() %>%
  bind_rows()

ELISpot_P_ST3 <-
sum_elispot  %>%
  as_tibble() %>%
  dplyr::select(contrast, estimate, p.value) %>%
  mutate(estimate = 10 ^ estimate) %>%
  mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value"
  ) %>%
  mutate(`Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
         `Fold change`= round(`Fold change`, digits = 6))%>%
 filter(`Adjusted P-value`< 0.05)

ELISpot_P_ST3 %>%
  DT::datatable(caption = "Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing interval") 
# Supplementary Table 3 visualization
t3v <-
  sum_elispot %>%
  as_tibble() %>%
  separate(contrast, c("a", "b"), sep = " - ") %>%
  mutate(a = str_remove_all(a , "interval"),
         b = str_remove_all(b , "interval")) %>%
  mutate(
    a = case_when(
      str_detect(a, "1") ~ paste(a, "week"),
      str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
      TRUE ~ paste(a, "weeks")
    ),
    b = case_when(
      str_detect(b, "1") ~ paste(b, "week"),
      str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
      TRUE ~ paste(b, "weeks")
    )
  ) %>%
  mutate(contrast = paste(a, "/", b)) %>%
  mutate(day = "1 week post-boost)") %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
  ) %>%
  mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
  ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
  geom_point(size = 3)  +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = .1,
                #position = position_dodge(width = 2),
                size = 1) +
  scale_color_manual(
    values = c(
      "#ABCFE5",
      "#93C4DE",
      "#78B5D8",
      "#60A6D1",
      "#4A97C9",
      "#3787C0",
      "#2575B7",
      "#1664AB",
      "#09539D",
      "#084185"
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
                     labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
  facet_grid( ~ day) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  labs(title = "S2-P-specific ASCs",
       y = "Estimated Fold-Change",
       x = "Prime-boost dosing interval") +
  coord_flip() +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    strip.text = element_text(size = 10),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

Model fitting for BM data

#-------------------------------------------------------------------------------
# Fit the model for BM data
#-------------------------------------------------------------------------------
# BM
d_bm <- subset(d, tissue == "BM")

boxcox_fits <-
  d_bm %>%
  split(.$coating) %>%
  map(~ boxcox_fit(.x))

lm_fits <-
  purrr::transpose(boxcox_fits)$lm_fit

# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
  emmeans(lm_grid, specs = ~ interval | day, type = "response") %>%
  as_tibble()

Residual diagnostics for BM data

assumption_plots <-
  imap(
    lm_fits ,
    ~ assumption_check(.x, .x$model$`linkfun(counts)`) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0("BM ", .y, "-specific model residual diagnostics "),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )
ggarrange(plotlist = assumption_plots)

Figure 4. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells

Each dot in the bottom panel of Figure 4 corresponds to the estimated mean for the prime-boost dosing interval. Error bars are 95% CIs on the estimated mean.

#-------------------------------------------------------------------------------
# Figure 4B. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
foldchange2 <-
  emmt_dat %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(
    factor(day),
    "4 week post-boost" = "87",
    "24 week post-boost" = "226"
  )) %>%
  mutate(day = factor(day, c("4 week post-boost" , "24 week post-boost"))) %>%
  ggplot(aes(x = interval , y = response, color = interval)) +
  geom_point(position = position_dodge(width = 0.5), size = 2)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 0.5),
    size = 1
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  labs(x = "", y = "Estimated counts/million cells") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(6)],
             color = "darkcyan",
             lty = "dashed") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(13)],
             color = "darkcyan",
             lty = "dashed") +
  facet_grid(~ day) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_blank(),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

box2 <-
  dp %>%
  filter(interval != "PBS", tissue == "BM", coating == "S2P") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(
    factor(day),
    "4 week post-boost" = "87",
    "24 week post-boost" = "226"
  )) %>%
  mutate(day = factor(day,
                      levels = c("4 week post-boost", "24 week post-boost"))) %>%
  ggplot(aes(x = interval, y = counts, color = interval)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "grey",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "Counts/million cells",
    color = "",
    linetype = "",
    title = "S2-P-specific LLPCs"
  ) +
  facet_grid(~ day) +
  theme(
    axis.text = element_text(size = 12),
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold")
  )

cowplot::plot_grid(
  box1,
  box2,
  foldchange1,
  foldchange2,
  labels = c("A", "C", "B", "D"),
  align = 'hv',
  rel_widths = c(1, 2),
  rel_heights = c(1, 1)
)

Supplementary Table 4. Statistical comparisons of S2-P–specific LLPCs 4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273 (10 μg) dosing intervals.

#-------------------------------------------------------------------------------
# Supplementary Table 4
#-------------------------------------------------------------------------------
emm_elispot_BM <-
  emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
          specs = ~ interval |
            day,
          type = "response")

set.seed(100)
emm_elispot_BM <-
  contrast(
    regrid(emm_elispot_BM , transform = "log10"),
    "pairwise",
    infer = TRUE,
    adjust = "mvt",
    ratio = FALSE,
    level = .95
  ) %>% as_tibble() %>%
  bind_rows()

ELISpot_P_ST4 <- 
emm_elispot_BM %>%
  as_tibble() %>%
  dplyr::select(day, contrast, estimate, p.value) %>%
  mutate(day = as.numeric(as.character(day))) %>%
  arrange(day) %>%
  mutate(estimate = 10 ^ estimate) %>%
  mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
    "Day" = "day"
  ) %>%
  mutate(`Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
         `Fold change`= round(`Fold change`, digits = 6))%>%
 filter(`Adjusted P-value`< 0.05)

ELISpot_P_ST4  %>%
  DT::datatable(caption = "Statistical comparisons of S2-P–specific LLPCs 4 weeks following
                dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273
                (10 μg) dosing intervals") 

Supplementary Figure 5. Statistical comparisons of S2-P-specific Antibody Secreting Cells and Long-Lived Plasma Cells between mRNA-1273 (10 µg) dosing intervals(Supplementary Table 3 and 4 visualization)

Each dot corresponds to the estimated fold change for the contrasted prime-boost dosing interval. Error bars are multivariate t adjusted 95% CIs on the estimated contrast.

t4v <-
  emm_elispot_BM %>%
  as_tibble() %>%
  separate(contrast, c("a", "b"), sep = " - ") %>%
  mutate(a = str_remove_all(a , "interval"),
         b = str_remove_all(b , "interval")) %>%
  mutate(
    a = case_when(
      str_detect(a, "1") ~ paste(a, "week"),
      str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
      TRUE ~ paste(a, "weeks")
    ),
    b = case_when(
      str_detect(b, "1") ~ paste(b, "week"),
      str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
      TRUE ~ paste(b, "weeks")
    )
  ) %>%
  mutate(contrast = paste(a, "/", b)) %>%
  mutate(day = factor(day, levels = c(87, 226))) %>%
  mutate(day = forcats::fct_recode(
    day,
    "4 week post-boost" = "87",
    "24 week post-boost" = "226",
  )) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
  ) %>%
  mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
  ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
  geom_point(size = 3)  +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = .1,
                #position = position_dodge(width = 2),
                size = 1) +
  scale_color_manual(
    values = c(
      "#ABCFE5",
      "#93C4DE",
      "#78B5D8",
      "#60A6D1",
      "#4A97C9",
      "#3787C0",
      "#2575B7",
      "#1664AB",
      "#09539D",
      "#084185"
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
                     labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  facet_grid(~ day) +
  labs(title = "S2-P-specific LLPC",
       # subtitle = "4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227)
       # between mRNA-1273 (10 μg) dosing intervals",
       y = "Estimated Fold-Change",
       x = "Prime-boost dosing interval") +
  coord_flip() +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    strip.text = element_text(size = 10),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

cowplot::plot_grid(
  plotlist = list(t3v, t4v),
  labels = c("A", "B"),
  nrow = 1,
  ncol = 2,
  rel_widths = c(0.5, 1)
)

Session information

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MASS_7.3-58.1   cowplot_1.1.1   ggpubr_0.4.0    DT_0.24        
##  [5] lme4_1.1-30     Matrix_1.5-1    ggeffects_1.1.3 emmeans_1.8.0  
##  [9] mgcv_1.8-40     nlme_3.1-159    readxl_1.4.1    forcats_0.5.2  
## [13] stringr_1.4.1   dplyr_1.0.10    purrr_0.3.4     readr_2.1.2    
## [17] tidyr_1.2.0     tibble_3.1.8    ggplot2_3.3.6   tidyverse_1.3.2
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.2            lubridate_1.8.0     httr_1.4.4         
##  [4] tools_4.2.1         backports_1.4.1     bslib_0.4.0        
##  [7] utf8_1.2.2          R6_2.5.1            DBI_1.1.3          
## [10] colorspace_2.0-3    withr_2.5.0         tidyselect_1.1.2   
## [13] compiler_4.2.1      cli_3.3.0           rvest_1.0.3        
## [16] xml2_1.3.3          labeling_0.4.2      sass_0.4.2         
## [19] scales_1.2.1        mvtnorm_1.1-3       digest_0.6.29      
## [22] minqa_1.2.4         rmarkdown_2.16      pkgconfig_2.0.3    
## [25] htmltools_0.5.3     highr_0.9           dbplyr_2.2.1       
## [28] fastmap_1.1.0       htmlwidgets_1.5.4   rlang_1.0.5        
## [31] rstudioapi_0.14     jquerylib_0.1.4     generics_0.1.3     
## [34] farver_2.1.1        jsonlite_1.8.0      crosstalk_1.2.0    
## [37] car_3.1-0           googlesheets4_1.0.1 magrittr_2.0.3     
## [40] Rcpp_1.0.9          munsell_0.5.0       fansi_1.0.3        
## [43] abind_1.4-5         lifecycle_1.0.1     stringi_1.7.8      
## [46] yaml_2.3.5          carData_3.0-5       grid_4.2.1         
## [49] crayon_1.5.1        lattice_0.20-45     haven_2.5.1        
## [52] splines_4.2.1       hms_1.1.2           knitr_1.40         
## [55] pillar_1.8.1        boot_1.3-28         estimability_1.4.1 
## [58] ggsignif_0.6.3      reprex_2.0.2        glue_1.6.2         
## [61] evaluate_0.16       modelr_0.1.9        vctrs_0.4.1        
## [64] nloptr_2.0.3        tzdb_0.3.0          cellranger_1.1.0   
## [67] gtable_0.3.1        assertthat_0.2.1    cachem_1.0.6       
## [70] xfun_0.32           xtable_1.8-4        broom_1.0.1        
## [73] coda_0.19-4         rstatix_0.7.0       googledrive_2.0.0  
## [76] gargle_1.2.0        writexl_1.4.0       ellipsis_0.3.2